function IRF = Gensys_IRF(G1, n_sectors, n_periods, pos_shock, imp)
%IRF: y(t+s) - E_t(y(t+s)) = sum_(v=0)^(s-1) G1^v*imp*z(t+s-v)

zt_imp = zeros(2 + n_sectors, n_periods);
zt_imp(pos_shock, 1) = 1;
IRF = zeros(size(G1, 2), n_periods + 1);   
for j = 1:n_periods
    IRF(:, j + 1) = G1*IRF(:,j) + imp*zt_imp(:,j);
end
